function [AV, AV_Alt, InitialShock_RelativeSize, InitialShock_RelativeSize_i, ModelReturnE_i, Loss_i, S, ModelReturn_j, LiquidationVolume_j,Fundlevel_Stress_i,LiquidationVolume_i, Loss_ShareCrossHoldings_i, ZeroPrice, TotalLiquidation, TotalLiquidation_NonFlow, TotalLiquidation_Flow] = AV_Calculation(TotalAssets_NON, TotalAssets_MKT, PortfolioMatrix_0, TotalAssets_FundShares, a_FundShares, Leverage, GammaEquity, Beta, Theta, PriceImpact, F_MKT, F_NON, Delta_log_VIX, fund_share_liquidations, onlyAV) % ,IV_FF, DV, utype

%%% This function computes the (aggregate) vulnerability from Fricke & Wilke (2023, RFS).

%%%%%%%%%%%%%%%%%%%

TotalAssets = TotalAssets_MKT + TotalAssets_NON + TotalAssets_FundShares;

%%% Compute debt and equity

Debt    = (Leverage ./ (1 + Leverage)) .* TotalAssets;
Equity  = TotalAssets - Debt;

%%% Step 1: Get matrix of investment weights. For each row, we now see the weight of each asset in that investor's portfolio

M = PortfolioMatrix_0./repmat(sum(PortfolioMatrix_0')',1,size(PortfolioMatrix_0,2));
M(isnan(M)) = 0;
a_FundShares = sparse(a_FundShares);

N = size(M, 1); % number of investors
K = size(M, 2); % number of marketable assets (including cash)

%%% Step 2: Get diagonal matrix of portfolio size (marketable assets, A_MKT, and all assets, A).

A       = sparse(diag(TotalAssets));
A_MKT   = sparse(diag(TotalAssets_MKT));
A_NON   = sparse(diag(TotalAssets_NON));
A_FUND  = sparse(diag(TotalAssets_FundShares));

%%% Check whether crossholdings exist
crossFundHoldingsExist = nnz(A_FUND) > 0;

%%% Step 3: Get diagonal matrix of target leverage (i. e., the pre-shock leverage), B.

B = sparse(diag(Leverage));

% Get diagonal matrix of equity
E     = sparse(diag(Equity));
E_agg = sum(Equity);

%%% Step 4: Get diagonal matrix of flow-sensitivities

Gamma = sparse(diag(GammaEquity));

%%% Step 5: Get diagonal matrix of price impacts

L = sparse(diag(PriceImpact));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% STEP 0 : INITIAL SHOCK %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Step 0: Initial Shock (direct effect)

%% Add shock = 0 for cash

if length(F_MKT) == 1 % for scalar input, generate shock vector (cash always the last asset)
    F_MKT = repmat(F_MKT, K, 1);
end

if PriceImpact(end) == 0
    F_MKT(end) = 0;
end

% Compute fund-level shocks to marketable and non-marketable assets
R0_MKT     = M * F_MKT; R0_MKT(isnan(R0_MKT)==1) = 0; % Eq. (1)
R0_NON     = F_NON;     % Eq. (2)

% Based on these, compute fund-level initial shocks
Btilde     = eye(N) + B;

R0_A       = (TotalAssets_MKT ./ TotalAssets) .* R0_MKT + (TotalAssets_NON ./ TotalAssets) .* R0_NON; % Eq. (6)
if crossFundHoldingsExist == 0
    R0_A(isnan(R0_A)) = 0;
end
R0_E       = max(diag(Btilde) .* R0_A, -1); % make sure the initial shock does not lead to negative equity.

InitialShock_RelativeSize_i = full(R0_E);
InitialShock_RelativeSize   = sum(full(R0_A .* diag(A))) / E_agg;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% STEP 1 : LossCrossHoldings_1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Step 1: Initial Shock (cross-holdings)

Katz_centrality = eye(N);
R1_E = R0_E;
R1_A = R0_A;
if crossFundHoldingsExist == 1
    Katz_centrality = inv(eye(N) - inv(E) * a_FundShares);
    R1_E       = Katz_centrality * R0_E;
    R1_E       = max(R1_E, -1);
    R1_A       = (1 ./ (1 + diag(B))) .* R1_E;
end

% End result: initial shock
R1A = R1_A;
R1E = R1_E;

InitialShockWithCrossHoldings_Euro_i = sum(full(R1_A .* diag(A)));
Loss_CrossHoldings_1                 = full(R1_A .* diag(A)) - full(R0_A .* diag(A));

% Compute Post-Shock-Portfolioholdings (PortfolioMatrix_1)
A_1               = A *(1 + R1A);
PortfolioMatrix_1 = PortfolioMatrix_0 .* repmat(1 + F_MKT', N, 1);
M_1               = PortfolioMatrix_1 ./ repmat(sum(PortfolioMatrix_1')', 1, K);
a_FundShares_1    = a_FundShares .* repmat(1 + R1E', N, 1);
M_0               = M;
M                 = M_1; % Replace initial portfolio weights matrix, M, by post-shock matrix, M_1.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%% STEP 2 : Flow-driven andn discrectionary (non-flow driven) liquidations %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Steps 2a and 2b

% Compute ...
... vector of fund-level overall amount to be liquidated, Phi
    Phi_Flow    = A * Gamma * R1A .* (1 + diag(Btilde).* R1A); % Eq. (21)
tmp_Flow    = Gamma*(diag(Btilde).*R1A);
Phi_NonFlow = A*(1 + R1A) .*(Beta.*tmp_Flow + Theta.*Delta_log_VIX);
Phi = Phi_Flow + Phi_NonFlow;

... vector of fund-level amount to be liquidated due to net outflows, Phi_NetOutflows
    Phi_Flow              = Phi - Phi_NonFlow;
Phi_ShareFlow         = abs(Phi_Flow ./ Phi);
Phi_ShareFlow(Phi==0) = 0;

Phi = max(full(-diag(A_MKT)), Phi); % Make sure that funds sell at most their marketable assets (not their total assets)

Phi_Flow     = Phi_ShareFlow .* Phi;
Phi_NonFlow  = (1 - Phi_ShareFlow) .* Phi;

LiquidationVolume_i = abs(Phi);
Fundlevel_Stress_i  = [abs(Phi) >= PortfolioMatrix_0(:, end) abs(Phi) PortfolioMatrix_0(:, end)];

... total amount to be liquidated by fund sector
    TotalLiquidation    = sum(Phi(:));

... total amount of discretionary liquidations
    TotalLiquidation_NonFlow = sum(Phi_NonFlow(:));

... total amount of flow-driven liquidations
    TotalLiquidation_Flow = sum(Phi_Flow(:));

%%%%%% now: distribute total amount in a pro-rata fashion

if fund_share_liquidations == 1
    M_tilde         = [PortfolioMatrix_1 a_FundShares_1];
    tmp_holdings    = sum(M_tilde')';
    M_tilde         = M_tilde ./ repmat(tmp_holdings, 1, N + K);
    M_tilde(isnan(M_tilde)) = 0;
    Z               = sparse([zeros(N, K) eye(N)]);
else
    M_tilde         = [PortfolioMatrix_1 zeros(N)];
    tmp_holdings    = sum(M_tilde')';
    M_tilde         = M_tilde ./ repmat(tmp_holdings, 1, N + K);

end

% Compute liquidation amount at the asset level
psi_init     =      M_tilde' * Phi;
psi_init_Flow =     M_tilde' * Phi_Flow;
psi_init_NonFlow  = M_tilde' * Phi_NonFlow;

if crossFundHoldingsExist == 1 & fund_share_liquidations == 1
    W       = sparse(eye(N + K) - M_tilde' * Z); % pre-allocate to speed up process
    psi     = W\psi_init;
    psi_Flow = W\psi_init_Flow;
    psi_NonFlow  = W\psi_init_NonFlow;
else
    psi = psi_init;
    psi_Flow= psi_init_Flow;
    psi_NonFlow = psi_init_NonFlow;
end

psi          = psi(1 : K);
psi_Flow     = psi_Flow(1 : K);
psi_NonFlow  = psi_NonFlow(1 : K);

psi_FundShareLiquidation  = psi - psi_init(1 : K);          % these are the bond + equity sales due to fund share redemptions
psi_AssetSales            = psi - psi_FundShareLiquidation; % these are the remaining bond + equity sales

LiquidationVolume_j = full(psi);
TotalLiquidation(2) = sum(Phi(:));
E2b = (diag(E) .* (1 + R1E)) + Phi + psi_init(K+1:end);
A2b = A_1 - (diag(E) - E2b);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% STEP 3 : Price Impact %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Step 3a: Direct losses due to price changes

% Compute �-price impact on assets sold by fire-selling funds
F3a_MKT             = max(L * psi, -1); % Eq. (14)
F3a_MKT_Flow        = F3a_MKT * abs(sum(Phi_Flow)) / abs(sum(Phi));
F3a_MKT_NonFlow     = F3a_MKT * abs(sum(Phi_NonFlow)) / abs(sum(Phi));

F3a_MKT_AssetSales           = F3a_MKT * sum(abs(psi_AssetSales)) / sum(abs(psi));
F3a_MKT_FundShareLiquidation = F3a_MKT * sum(abs(psi_FundShareLiquidation)) / sum(abs(psi));

ModelReturn_j = [full(F3a_MKT) full(F3a_MKT_AssetSales) full(F3a_MKT_FundShareLiquidation)]; % asset level returns

... finally, compute %-price impact on assets sold by fire-selling funds
    R3a_MKT          = M * F3a_MKT;
R3a_MKT_Flow     = M * F3a_MKT_Flow;
R3a_MKT_NonFlow  = M * F3a_MKT_NonFlow;
R3a_MKT_AssetSales              = M * F3a_MKT_AssetSales;
R3a_MKT_FundShareLiquidation    = M * F3a_MKT_FundShareLiquidation;

% Direct losses
Loss_FireSales = diag(A_MKT) .* R3a_MKT;
Loss_FireSales = full(max(-diag(A), Loss_FireSales)); % funds cannot lose more than their total assets

w_Flow                      = nansum(diag(A_MKT) .* R3a_MKT_Flow) / (nansum(diag(A_MKT) .* R3a_MKT_Flow) + nansum(diag(A_MKT) .* R3a_MKT_NonFlow));
Loss_FireSales_Flow         = full(w_Flow * Loss_FireSales);
Loss_FireSales_NonFlow      = full((1 - w_Flow) * Loss_FireSales);

w_AssetSales                = nansum(diag(A_MKT).* R3a_MKT_AssetSales) / (nansum(diag(A_MKT) .* R3a_MKT_AssetSales) + nansum(diag(A_MKT) .* R3a_MKT_FundShareLiquidation));
if crossFundHoldingsExist == 0
    w_AssetSales = 1;
end
Loss_FireSales_AssetSales           = full(w_AssetSales * Loss_FireSales);
Loss_FireSales_FundShareLiquidation = full((1 - w_AssetSales) * Loss_FireSales);

R3a_E                                   = (Loss_FireSales) ./ E2b;
R3a_E(isnan(R3a_E) & R1_E == -1)        = 0;
R3a_E(isinf(abs(R3a_E)) & R1_E == -1)   = 0;

if crossFundHoldingsExist == 0
    R3a_E(isnan(R3a_E)) = 0;
end

%%% Step 3b: Cross-holdings

R3b_E                = Katz_centrality * R3a_E;
Katz_centrality      = (Katz_centrality * ones(N, 1) - 1); % normalize Katz centrality such that funds with no indegree have zero Katz
if crossFundHoldingsExist == 0
    Katz_centrality = Katz_centrality + 1;
end

Loss_CrossHoldings_3 = (R3b_E - R3a_E) .* diag(full(E)); 
Loss_CrossHoldings   = Loss_CrossHoldings_1 + Loss_CrossHoldings_3;

ModelReturnE_i            = (Loss_FireSales + Loss_CrossHoldings) ./ diag(E);
Loss_ShareCrossHoldings_i = Loss_CrossHoldings ./ (Loss_FireSales + Loss_CrossHoldings);
if crossFundHoldingsExist == 0
    ModelReturnE_i(isnan(ModelReturnE_i)) = 0;
    Loss_ShareCrossHoldings_i = zeros(length(Loss_ShareCrossHoldings_i),1);
end

%%%%%%%%%%%%%% Aggregate Vulnerability (AV)

% Compute the fund sector's overall vulnerability, relative to its total pre-shock size

AV_numerator_Flow            = sum(Loss_FireSales_Flow);
AV_numerator_NonFlow         = sum(Loss_FireSales_NonFlow);
AV_numerator_CrossHoldings   = sum(Loss_CrossHoldings);

% Baseline aggregate vulnerability

AV                       = (AV_numerator_Flow + AV_numerator_NonFlow + AV_numerator_CrossHoldings) / E_agg;
AV_Flow                  = AV_numerator_Flow / E_agg;
AV_NonFlow               = AV_numerator_NonFlow / E_agg;
AV_CrossHoldings         = AV_numerator_CrossHoldings / E_agg;
AV_CrossHoldings_Initial = (AV_numerator_CrossHoldings - sum(Loss_CrossHoldings_3)) / E_agg;

AV = [AV AV_Flow AV_NonFlow AV_CrossHoldings AV_CrossHoldings_Initial];

AV_numerator_CrossHoldings_1                 = sum(Loss_CrossHoldings_1);
AV_numerator_CrossHoldings_3                 = sum(Loss_CrossHoldings_3);
AV_numerator_FireSales_AssetSales            = sum(Loss_FireSales_AssetSales);
AV_numerator_FireSales_FundShareLiquidation  = sum(Loss_FireSales_FundShareLiquidation);

AV_CrossHoldings_1                  = AV_numerator_CrossHoldings_1 / E_agg;
AV_FireSales_AssetSales             = AV_numerator_FireSales_AssetSales / E_agg;
AV_FireSales_FundShareLiquidation   = AV_numerator_FireSales_FundShareLiquidation / E_agg;
AV_CrossHoldings_3                  = AV_numerator_CrossHoldings_3 / E_agg;

AV_Alt = [AV(1) AV_FireSales_AssetSales AV_FireSales_FundShareLiquidation AV_CrossHoldings_1 AV_CrossHoldings_3];

ZeroPrice                = zeros(K, 1);
ZeroPrice(F3a_MKT == -1) = 1;

%%% Individual values

S             = NaN(N, 3);
Loss_i        = [Loss_CrossHoldings_1 Loss_FireSales_AssetSales Loss_FireSales_FundShareLiquidation Loss_CrossHoldings_3]./repmat(diag(E),1,4);

% Compute each fund's systemicness (S)
if onlyAV == 0
    S_SalesDirect = zeros(N, 1);
    S_CrossFund   = (Katz_centrality ./ sum(Katz_centrality)) .* sum(Loss_CrossHoldings);% katz centrality
    for i = 1:N
        % Here, we attribute the total price changes in Step 3a proportional to the asset sales of fund f in each asset.
        phi_i = M_tilde(i,:)*Phi(i);
        phi_i = phi_i(1:K)';
        phi   = psi_AssetSales(1:K);
        sh_i  = phi_i./phi;
        sh_i(phi == 0) = 0;

        S_SalesDirect(i)  = sum(PortfolioMatrix_1*(F3a_MKT_AssetSales.*sh_i));
    end

    % Normalize systemicness (if necessary)
    if sum(Loss_CrossHoldings) < 0
        S_CrossFund = S_CrossFund./E_agg;
        S_CrossFund = sum(AV_Alt(3:end))/sum(S_CrossFund).*S_CrossFund;
    end
    S_SalesDirect = S_SalesDirect./E_agg;
    S_SalesDirect = AV_Alt(2)/sum(S_SalesDirect).*S_SalesDirect;

    S = S_CrossFund + S_SalesDirect;
    S = [S S_CrossFund S_SalesDirect];
end